      SUBROUTINE CTPAK3A(KSYSTEM,KORIGIN,ROTMTX,ORIGIN50,
     *             NPARMSIN,PARMSIN,LUERR,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C PURPOSE:  COMPUTE MATRIX TO ROTATE VECTORS FROM A SPECIFIED SYSTEM TO
C           MEAN EARTH EQUATOR AND EQUINOX OF 1950.0, AND COMPUTE 
C           TRANSLATION VECTOR TO SHIFT THE ORIGIN TO GEOCENTRIC, MEAN 
C           OF 1950.0
C
C***********************************************************************
C
C  BY K PACKARD, 4/83.
C      MODIFICATIONS - EXTENSIVE MODS BY C PETRUZZO, 10/83, 7/84.
C
C***********************************************************************
C
C   KSYSTEM, ASSIGNMENTS FOR THE 'FROM' SYSTEM'S ORIENTATION :
C
C     1  MEAN EARTH EQUATOR AND EQUINOX OF 1950.0
C     2  MEAN EARTH EQUATOR AND EQUINOX OF DATE
C     3  TRUE EARTH EQUATOR AND EQUINOX OF DATE
C     4  MEAN EARTH EQUATOR AND EQUINOX OF 2000.0
C
C     5  MEAN ECLIPTIC AND EQUINOX OF 1950.0
C     6  MEAN ECLIPTIC AND EQUINOX OF DATE
C     7  TRUE ECLIPTIC AND EQUINOX OF DATE. NOT AVAILABLE. DOES MEAN.
C     8  MEAN ECLIPTIC AND EQUINOX OF 2000.0
C
C     9  GALACTIC. INT'NL ASTRO UNION CONVENTION, 1958.
C    10  TRUE EARTH EQUATOR, X-AXIS THROUGH GREENWICH MERIDIAN
C    11  TOPOCENTRIC, X-AXIS SOUTH
C    12  LOCAL ORBITAL
C    13  SPACECRAFT BODY
C    14  ORBIT PLANE, X-AXIS THROUGH ASCENDING NODE
C
C   KORIGIN, ASSIGNMENTS FOR THE 'FROM' SYSTEM'S ORIGIN :
C
C        1 -> EARTH CENTER
C        2 -> SUN CENTER, MEAN OF 1950.0
C        3 -> SPACECRAFT CENTER, MEAN OF 1950.0
C        4 -> TOPOCENTRIC
C
C
C
      REAL*8 HALFPI/ 1.570796326794897D0 /
      REAL*8 PI/ 3.141592653589793D0 /
      REAL*8 TWOPI/ 6.283185307179586D0 /
      REAL*8 DEGRAD/ 57.29577951308232D0 /
C
      REAL*8 ROTMTX(3,3),SCPOS(3),SCVEL(3),ORIGIN50(3)
      REAL*8 ORIGINTOP(3),ROT1(3,3),ROT2(3,3),ROTGEOGR(3,3)
      REAL*8 IDENT(3,3)/1.D0,3*0.D0,1.D0,3*0.D0,1.D0/
      REAL*8 ROT2000(3,3),ROTECL2000(3,3),ROTECL1950(3,3),ROTGALAC(3,3)
      LOGICAL NEEDGEOGR,NEED2000
C
      INTEGER INIT2000/0/,INIT5/0/,INIT8/0/,INIT9/0/
C
      PARAMETER NPARMS = 16    ! MUST EQUAL NPARMSIN
      REAL*8 PARMS(NPARMS),PARMSIN(NPARMSIN)
      EQUIVALENCE            (PARMS(1),RATTREF),
     *   (PARMS(2),ZRA),     (PARMS(3),ZDEC),     (PARMS(4),CLOCK),
     *   (PARMS(5),TIME),    (PARMS(6),SCPOS(1)), (PARMS(9),SCVEL(1)), 
     *   (PARMS(12),TOPLAT), (PARMS(13),TOPLONG), (PARMS(14),TOPALT),
     *   (PARMS(15),FLAT),   (PARMS(16),ERAD)
C
C
C
      IF(NPARMS.NE.NPARMSIN) STOP ' STOPPED IN CTPAK3A. SOURCE ERROR.'
C
C
C
      IBUG=0
      LUBUG=19
C
C
C INITIALIZATION FOR THIS CALL.
      IERR = 0
      CALL MTXEQL(PARMSIN,PARMS,NPARMS,1)   ! SEE EQUIVALENCE STMT ABOVE
      IATTREF=JIDNNT(RATTREF)
C
C
C DEBUG
      IF(IBUG.NE.0) WRITE(LUBUG,9549) PARMS
 9549 FORMAT(/,' CTPAK3A. PARMS='/,4(4G18.9/))
      IF(IBUG.NE.0) WRITE(LUBUG,9550) KSYSTEM,KORIGIN,
     *     IATTREF,ZRA*DEGRAD,ZDEC*DEGRAD,CLOCK*DEGRAD,PAKTIM50(TIME),
     *     SCPOS,SCVEL,TOPLAT*DEGRAD,TOPRACN*DEGRAD,TOPALT,FLAT,ERAD
 9550 FORMAT(/,' CTPAK3A. '/,
     *         '    KSYSTEM=',I2,'  KORIGIN=',I2/,
     *         '    IATTREF,ZRA,ZDEC,CLOCK= ',I5,3F11.2/,
     *         '    TIME(PACKED)=',G22.13/,
     *         '    SCPOS=',3G15.7/,
     *         '    SCVEL=',3G15.7/,
     *         '    TOPLAT,TOPRACN,TOPALT=',3G15.7/,
     *         '    FLAT=',G13.5,'  ERAD=',G14.7)
C
C
C *************************************
C * CODE FOR MULTIPLE NEED SITUATIONS *    AVOIDS REPETITION OF CODE
C *************************************
C
C    CHECK NEED TO HAVE A ROTATION MATRIX TO GO FROM GEOGRAPHIC
C    ORIENTATION TO MEAN OF 1950.0
C
      NEEDGEOGR = KSYSTEM.EQ.10  .OR.  KSYSTEM.EQ.11  .OR.
     *            KORIGIN.EQ.4
C
      IF(NEEDGEOGR) THEN
        GHA = GHANGL(TIME,1)   !MEAN OF DATE GHA; T.O.D. WOULD BE BETTER
C      GET ROT1, FROM GEOGRAPHIC TO T.O.D.
        CALL ROT123(3,-GHA,0,DUM,0,DUM,ROT1,0,KDUM)
C      GET ROT2, FROM T.O.D. TO MEAN OF 1950.0
        CALL M50TOD(TIME,-1,ROT2)
C      FINALLY, FROM GEOGRAPHIC TO MEAN OF 1950.0
        CALL MTXMUL33(1,ROT2,ROT1,ROTGEOGR)
        END IF
C
C    CHECK NEED TO HAVE A ROTATION MATRIX TO GO FROM MEAN OF 2000
C    ORIENTATION TO MEAN OF 1950.0
C
      NEED2000 = KSYSTEM.EQ.4 .OR. KSYSTEM.EQ.8
C
      IF(NEED2000) THEN
        IF(INIT2000.EQ.0) THEN
          INIT2000 = 1
          TEMP = SINCE50(991231.235959D0) + 1.D0
          CALL M50MDT(TEMP,-1,ROT2000)
          END IF
        END IF
C
C
C
C *******************************
C * ROTATION MATRIX COMPUTATION *
C *******************************
C
C
      GO TO (  100,  200,  300,  400,  500,
     *         600,  700,  800,  900, 1000,
     *        1100, 1200, 1300, 1400),      KSYSTEM
C
C
C  >>> FROM MEAN EQUATOR AND EQUINOX OF 1950.0 <<<
C
  100 CONTINUE
      CALL MTXEQL(IDENT,ROTMTX,3,3)
      GO TO 5000
C
C
C  >>> FROM MEAN EQUATOR AND EQUINOX OF DATE <<<
C
  200 CONTINUE
      CALL M50MDT(TIME,-1,ROTMTX)
      GO TO 5000
C
C
C  >>> FROM TRUE EQUATOR AND EQUINOX OF DATE <<<
C
  300 CONTINUE
      CALL M50TOD(TIME,-1,ROTMTX)
      GO TO 5000
C
C
C  >>> FROM MEAN EQUATOR AND EQUINOX OF 2000.0 <<<
C
  400 CONTINUE
      CALL MTXEQL(ROT2000,ROTMTX,3,3)
      GO TO 5000
C
C
C  >>> FROM MEAN ECLIPTIC AND EQUINOX OF 1950.0 <<<
C
  500 CONTINUE
      IF(INIT5.EQ.0) THEN
        INIT5 = 1
        TEMP = OBLIQ(0.D0)
        CALL ROT123(1,-TEMP,0,DUM,0,DUM,ROTECL1950,0,KDUM)
        END IF
      CALL MTXEQL(ROTECL1950,ROTMTX,3,3)
      GO TO 5000
C
C
C  >>> FROM MEAN ECLIPTIC AND EQUINOX OF DATE <<<
C
  600 CONTINUE
      CALL ROT123(1,-OBLIQ(TIME),0,DUM,0,DUM,ROT1,0,KDUM)
      CALL M50MDT(TIME,-1,ROT2)
      CALL MTXMUL33(1,ROT2,ROT1,ROTMTX)
      GO TO 5000
C
C
C  >>> FROM TRUE ECLIPTIC AND EQUINOX OF DATE <<<
C
  700 CONTINUE
      GO TO 600                         ! <<--------------------
C
C
C  >>> FROM MEAN ECLIPTIC AND EQUINOX OF 2000.0 <<<
C
  800 CONTINUE
      IF(INIT8.EQ.0) THEN
        INIT8 = 1
        TEMP = SINCE50(991231.235959D0) + 1.D0
        CALL ROT123(1,-OBLIQ(TEMP),0,DUM,0,DUM,ROT1,0,KDUM)
        CALL MTXMUL33(1,ROT2000,ROT1,ROTECL2000)
        END IF
      CALL MTXEQL(ROTECL2000,ROTMTX,3,3)
      GO TO 5000
C
C
C  >>> FROM GALACTIC <<<
C
  900 CONTINUE
      IF(INIT9.EQ.0) THEN
        INIT9 = 1
        CALL VECM50GAL(-1,IDENT(1,1),ROTGALAC(1,1))
        CALL VECM50GAL(-1,IDENT(1,2),ROTGALAC(1,2))
        CALL VECM50GAL(-1,IDENT(1,3),ROTGALAC(1,3))
        END IF
      CALL MTXEQL(ROTGALAC,ROTMTX,3,3)
      GO TO 5000
C
C
C  >>> FROM GEOGRAPHIC <<<
C
 1000 CONTINUE
      CALL MTXEQL(ROTGEOGR,ROTMTX,3,3)
      GO TO 5000
C
C
C  >>> FROM TOPOCENTRIC <<<
C
 1100 CONTINUE
C    ROTATE FROM TOPOCENTRIC TO GEOGRAPHIC
      YAW = TOPLONG
      PITCH = HALFPI - TOPLAT
      CALL ROT123(3,YAW,2,PITCH,0,DUM,ROT1,0,KDUM)
C    ROTATE FROM GEOGRAPHIC TO MEAN OF 1950.0
      CALL MTXMUL33(2,ROTGEOGR,ROT1,ROTMTX)
      GO TO 5000
C
C
C  >>> FROM LOCAL ORBITAL <<<
C
 1200 CONTINUE
      CALL TRANS(SCPOS,SCVEL,ROTMTX)
      CALL MTXFLIP33(ROTMTX,ROTMTX)
      GO TO 5000
C
C
C  >>> FROM SPACECRAFT BODY <<<
C
 1300 CONTINUE
      CALL ATTMX2(-1,ZRA,ZDEC,CLOCK,IATTREF,ROTMTX,SCPOS,SCVEL)
      GO TO 5000
C
C
C  >>> FROM ORBIT PLANE, X-AXIS THROUGH ASCENDING NODE <<<
C
 1400 CONTINUE
      CALL UCROSS(SCPOS,SCVEL,ORBNORM)
      CALL XYZSPH(-1,AZNORM,ELNORM,DUM,ORBNORM,1.D0)
      ORBINCL = HALFPI - ELNORM
      ORBNODE = EQVANG(AZNORM+HALFPI)
      YAW = ORBNODE
      ROLL = ORBINCL
      CALL ROT123(3,YAW,1,ROLL,0,DUM,ROTMTX,0,KDUM)
      CALL MTXFLIP33(ROTMTX,ROTMTX)
      GO TO 5000
C
C
C
 5000 CONTINUE
C
C **********************
C * TRANSLATION VECTOR *
C **********************
C
C FIND THE COORDINATES OF THE 'FROM' SYSTEM'S ORIGIN IN THE
C MEAN OF 1950.0 SYSTEM.
C
C
      IF(KORIGIN.EQ.1) THEN
        ORIGIN50(1)=0.D0
        ORIGIN50(2)=0.D0
        ORIGIN50(3)=0.D0
        GO TO 5500
        END IF
C
      IF(KORIGIN.EQ.2) THEN 
        TIME=PARMS(5)
        CALL SOLM50(TIME,ORIGIN50,0)
        GO TO 5500
        END IF
C
      IF(KORIGIN.EQ.3) THEN
        ORIGIN50(1) = SCPOS(1)
        ORIGIN50(2) = SCPOS(2)
        ORIGIN50(3) = SCPOS(3)
        GO TO 5500
        END IF
C
      IF(KORIGIN.EQ.4) THEN
C      CONVERT ORIGIN TO EARTH-FIXED, GEOGRAPHIC.
        CALL TOPCENXYZ(TOPLAT,TOPLONG,TOPALT,ERAD,FLAT,ORIGINTOP)
C      NOW CONVERT TO MEAN OF 1950.0
        CALL VECROT(ROTGEOGR,ORIGINTOP,ORIGIN50)
        GO TO 5500
        END IF
C
C
 5500 CONTINUE
C
C
C
      RETURN
      END
